
#include "kubatko_prism.hpp"

namespace neon
{
kubatko_prism::kubatko_prism(int const minimum_degree)
{
    switch (minimum_degree)
    {
        case 1:
        {
            m_degree = 1;
            m_weights = {4.0};
            m_coordinates = {{0, -1.0 / 3.0, -1.0 / 3.0, 0.0}};

            break;
        }
        case 2:
        {
            m_degree = 2;

            m_weights.resize(4, 1.0);

            m_coordinates = {{0, 0.483163247594393, -0.741581623797196, 0.0},
                             {1, -0.605498860309242, 0.469416096821288, 0.0},
                             {2, -0.605498860309242, -0.530583903178712, 0.816496580927726},
                             {3, -0.605498860309242, -0.530583903178712, -0.816496580927726}};

            break;
        }
        case 3:
        {
            m_degree = 3;

            m_weights = {0.742534890852309,
                         0.375143463443327,
                         0.495419047908462,
                         0.523999970843238,
                         0.980905839025611,
                         0.881996787927053};

            m_coordinates = {{0, 0.240692796349625, -0.771991660873412, 0.614747128207527},
                             {1, -0.968326281451138, -0.568046512457875, 0.676689529541421},
                             {2, 0.467917833640195, -0.549342790168347, -0.599905857322635},
                             {3, -0.786144119530819, 0.362655041695561, -0.677609795694786},
                             {4, -0.484844897886675, -0.707931130717342, -0.502482717716373},
                             {5, -0.559053711932125, 0.260243325246813, 0.493010512161538}};

            break;
        }
        case 4:
        {
            m_degree = 4;
            m_weights = {0.545658450421913,
                         0.545658450421913,
                         0.545658450421913,
                         0.431647899262139,
                         0.431647899262139,
                         0.249954808368331,
                         0.249954808368331,
                         0.249954808368331,
                         0.249954808368331,
                         0.249954808368331,
                         0.249954808368331};

            m_coordinates = {{0, -0.062688380276010, -0.062688380276010, 0.0},
                             {1, -0.062688380276010, -0.874623239447980, 0.0},
                             {2, -0.874623239447980, -0.062688380276010, 0.0},
                             {3, -1.0 / 3.0, -1.0 / 3.0, 0.866861974009030},
                             {4, -1.0 / 3.0, -1.0 / 3.0, -0.866861974009030},
                             {5, -0.798519188402179, -0.798519188402179, 0.675639823682265},
                             {6, -0.798519188402179, -0.798519188402179, -0.675639823682265},
                             {7, -0.798519188402179, 0.597038376804357, 0.675639823682265},
                             {8, -0.798519188402179, 0.597038376804357, -0.675639823682265},
                             {9, 0.597038376804357, -0.798519188402179, 0.675639823682265},
                             {10, 0.597038376804357, -0.798519188402179, -0.675639823682265}};
            break;
        }
        case 5:
        {
            m_degree = 5;

            m_weights = {0.213895020288765,
                         0.141917375616806,
                         0.295568859378071,
                         0.256991945593379,
                         0.122121979248381,
                         0.175194917962627,
                         0.284969106392719,
                         0.323336131783334,
                         0.159056110329063,
                         0.748067388709644,
                         0.280551223607231,
                         0.147734016552639,
                         0.259874920383688,
                         0.235144061421191,
                         0.355576942732463};

            m_coordinates = {{0, -0.820754415297359, 0.701020947925133, -0.300763696502910},
                             {1, 0.611831616907812, -0.869995576950389, 0.348546607420888},
                             {2, -0.951379065092975, -0.087091980055873, 0.150656042323906},
                             {3, 0.200535109198601, -0.783721735474016, -0.844285153176719},
                             {4, -0.909622841605196, -0.890218158063352, 0.477120081549168},
                             {5, 0.411514133287729, -0.725374126531787, 0.864653509536562},
                             {6, -0.127534496411879, -0.953467283619037, 0.216019762875977},
                             {7, -0.555217727817199, -0.530472194607007, 0.873409672725819},
                             {8, 0.706942532529193, -0.782248553944540, -0.390653804976705},
                             {9, -0.278092963133809, -0.291936530517119, -0.126030507204870},
                             {10, -0.057824844208300, -0.056757587543798, 0.539907869785112},
                             {11, -0.043308436222116, 0.012890722780611, -0.776314479909204},
                             {12, -0.774478920734726, 0.476188541042454, 0.745875967497062},
                             {13, -0.765638443571458, 0.177195164202219, -0.888355356215127},
                             {14, -0.732830649614460, -0.737447982744191, -0.651653242952189}};
            break;
        }
        case 6:
        {
            throw std::domain_error("Sixth order scheme not yet implemented");
            break;
        }
        default:
        {
            throw std::domain_error("This is not the quadrature scheme you are looking for");
        }
    }
    // Transform from their reference prism to ours

    // mapping for weights is (weight / 4.0)
    std::transform(begin(m_weights), end(m_weights), begin(m_weights), [](auto const weight) {
        return weight / 4.0;
    });

    // mapping for r and s is (xi + 1.0) / 2.0
    // mapping for z is (zeta + 1.0) / 2.0
    std::transform(begin(m_coordinates),
                   end(m_coordinates),
                   begin(m_coordinates),
                   [](auto coordinate_pack) {
                       auto& [i, x, y, z] = coordinate_pack;

                       x = (x + 1.0) / 2.0;
                       y = (y + 1.0) / 2.0;
                       z = (z + 1.0) / 2.0;

                       return coordinate_pack;
                   });
}
}
